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1. Introduction 



The observation of Bose-Einstein condensation in 
atomic traps [l]-[3] has renewed theoretical interest in 
studying a system of Bose particles with a nonuniform 
density. The experiments involve atoms of mass m con- 
densed in a harmonic trap characterized by frequency 
o). The relevant parameters, as discussed more fully in 
Refs. [4] and [5], have the following orders of magni- 
tude: 

a - 10"^ cm, 

A = {ln/i^lmkTY^ - 10"^ cm, 

L={hlmoyf'^- 10"^ cm. 



where a is the ^-wave scattering length, A the thermal 
wavelength, L the size of the system, and n is the average 
particle density. Thus, na^ « l,a/L« l,a/\« 1, and 
we can treat the system as a low- temperature weakly-in- 
teracting dilute gas using well-known methods [6]. 

In this paper, we review the following theoretical 
methods, and point out some common errors and mis- 
conceptions: 

(a) The pseudopotential method: One replaces the 
potential by a 5-function, and, with the help of a Bogoli- 
ubov transformation, obtains a weak-coupling expan- 



10" 



(1) 
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sion. This method seems adequate to describe current 
experiments. 

(b) The self-consistent field method: This is varia- 
tional, and capable of treating strongly interacting cases 
in principle. It reduces to the pseudopotential method in 
appropriate limits; but one should not use a 5 -function 
potential in a variational calculation, because such a 
potential in three spatial dimensions has no effect, if 
treated exactly. 

(c) The Gaussian variational method: One uses a 
Gaussian trial wave function, and obtains results similar 
to the self-consistent field method. It has the advantage 
of easily adaptable to the description of time-dependent 
phenomena. 

All three methods are closely related to one another. 
Although we will not use Feynman graphs, it might be 
relevant to note that the Bogoliubov transformation cor- 
responds to summing ''one-loop" graphs. The self-con- 
sistent field method consists of readjusting parameters 
in the one-loop sum in a variational sense. The Gaussian 
method includes all one-loop contribution, plus parts of 
higher-loop contributions. Of these methods, the 
straightforward one-loop summation yields an expan- 
sion in na^. On the other hand, the accuracy of the other 
methods cannot be ascertained. 

We do not attempt a comprehensive review, and the 
references cited are by no means complete. In a subject 
with such a long history, it is difficult to be aware of all 
relevant sources, and we apologize for any omissions [7]. 



2. Pseudopotential Method 

In low- temperature calculations, one can replace the 
interatomic potential with a 5 -function potential 
(47Ta^Vm)5(r), where a > is the ^-wave scattering 
length. Such a potential should be treated in low-order 
perturbation theory only, because the 5 -function poten- 
tial in three spatial dimensions has no effect if treated 
exactly [8]. 

Consider the Schrodinger equation 



-V'iP(r)^g8\r)iP(r) = EiP(rl 



(2) 



where E is 2m Iff' times the energy. Taking the Fourier 
transform of both sides leads to 



{E-e)4>{k) = gifj{Q\ 



(3) 



where (^(A:) is the Fourier transform of i/f(r). This must 
hold for all A:, including k^ = E. Therefore ij/(0) = 0. 
This condition is automatically satisfied for all partial 
waves except ^ -waves. For ^ -waves, this seems to re- 
quire that ij/(r) be discontinuous at r = 0; but because of 



the singular nature of the potential at r = 0, one should 
be careful. By treating the 5-function as the limit of a 
square well, one can verify that in three dimensions the 
^ -waves are the same as those of a free particle except 
at r = 0, where they drop to zero discontinuously. This 
behavior, however, does not lead to any scattering, nor 
level shift. 

As a more careful analysis [9] shows, the 5-function 
potential must be regarded as the first term in a pseudo- 
potential, which correctly reproduces all scattering 
phase shifts, and depends on these phase shifts as input 
parameters. For the simple case of a hard- sphere poten- 
tial, for which there is only one parameter a, the hard- 
sphere diameter, the pseudopotential has the form 



' pseudt 



4TTa;5^ 



,o(r): 



8^r)lr-^-8W^lr. 



(4) 



The differential operator {dldr)r removes spurious sin- 
gularities in the wave function, which would otherwise 
make the ground- state energy diverge. In low-order per- 
turbation theory, one can replace this operator by a sim- 
ple subtraction procedure for the energy. 

Using the first term in the pseudopotential, and mak- 
ing a Bogoliubov transformation, one can calculate ex- 
actly the first few terms of an expansion for the ground- 
state energy per particle of a uniform hard-sphere Bose 
gas [10, 11]: 



^0 = 



2iina^ 



m 






(f-V3) 



na^\rv{na^) + 0(na^) 



(5) 



This shows that the energy is not analytic at a = 0, and 
therefore cannot be continued to negative values of a . 

Fetter [12] has given a systematic treatment of the 
nonuniform case based on the Bogoliubov transforma- 
tion. For application to the atomic-trap experiments, one 
can also adapt the results of the uniform case using a 
Thomas-Fermi approximation [5]. 

What happens when the scattering length is negative? 
In this case, the two-body scattering is dominated by the 
attractive part of the potential. The attraction may or 
may not be strong enough to produce a two-particle 
bound state; but it will lead to an A/^-body bound state, 
for a sufficiently large N. The many-body problem is 
very different from that of the repulsive case, for it 
depends on the details of the potential. Consider two 
potentials with the same negative scattering length, one 
everywhere attractive, and the other having a repulsive 
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core. The many-body system with the purely attractive 
potential will collapse spatially, whereas the one with 
hard core will yield extensive energies. 

When the scattering length is negative, we must there- 
fore use an actual potential, for example, one with a hard 
core plus an attractive tail. We can replace the hard-core 
part by a pseudopotential, and do the Bogoliubov trans- 
formation. In this manner, it has been shown [6] [13] 
that a dilute Bose system with such a potential exhibits 
a first-order phase transition separating a gas phase 
from a liquid phase. The transition line terminates in a 
critical point. The Bose-Einstein condensation occurs in 
the liquid phase, and divides it into liquid I and liquid II. 
In short, one reproduces qualitatively the phase diagram 
of "^He. Since the calculation is perturbative, one can 
follow the collapse of the phase diagram to that of the 
ideal Bose gas when the potential is turned off. 

The results described above pertains to an infinite 
Bose system in the thermodynamic limit, and may be 
modified for a finite system. The authors of Ref. [2] 
observed Bose condensation in an atomic trap filled 
with in ^Li, which has a negative scattering length. In 
this case, the smallness of the system may have rendered 
the first-order transition indistinct. It is also possible that 
the observed phenomenon is metastable. 



3. Variational Methods and the Energy 
Gap 

A variational approach to the ground- state energy of 
a Bose gas was introduced by Girardeau and Arnowitt 
[14], who use (wisely) an actual potential instead of the 
5-function. Their results are for a uniform gas, but oth- 
erwise very similar to what we obtain later in the self- 
consistent field method. The excitation spectrum in this 
approach, however, contains an energy gap, contrary to 
general expectations [15], and to the perturbation-the- 
ory results in the hard-sphere case [10]. 

A variational method is designed to yield the best 
ground-state energy, and may not yield the excitation 
spectrum accurately. In this instance, we know that the 
gap should have been filled with phonons, which can be 
described using general principles. We may therefore 
view the variational principle as a way to obtain an 
approximate ground-state wave function, and build the 
phonon states upon this ground state. 

Feynman [16] argues that the low- lying excited states 
of a Bose system are purely phonon states, owing to the 
statistics, and suggests the following one-phonon wave 
function of momentum k : 



where ^o is the ground- state wave function. The excita- 
tion energy is 



Ek 



hV 



2mSk ' 
where St is the Uquid structure factor 



(7) 



(8) 



where 



i, = n-"'^a;^. 



(9) 



with ak the annihilation operator of a free particle of 
momentum M , and O the volume. From general princi- 
ples [17,18], we expect 



Sk-^ 



& 



^k^o 2mc 



(10) 



where c is the sound velocity as computed from the 
compressibility of the ground state. This leads to the 
phonon spectrum 



Ek = hck. 



(11) 



Feynman's argument can be sharpened by recasting it as 
the statement that the longitudinal sum-rules are satu- 
rated by the phonon states [17] [18]. All the equations 
above are verified in the dilute hard-sphere Bose gas 
[10] [18]. 

Other types of states that can be built upon ^o are 
those describing superfluid flow [10] [18]: 

^supcmow(ri,...,r^) = ne'-"^0-^^o(ri,.../^), (12) 

j 

where a{r) is the superfluid phase, which is related to 
the superfluid velocity v^ through 



Vs(r)=- Va(r). 
m 



(13) 



With this wave function one can describe quantized vor- 
tices. 



4. Self-Consistent Field Method: Ground 
State 



A^ 






'i^o{ru-,rN), 



(6) 



The self-consistent field method was first used by 
Bogoliubov in his treatment of superconductivity, and is 
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therefore also known as the "Bogoliubov method." We 
shall follow a concise formulation due to De Gennes 
[19]. A recent review of similar methods is given by 
Griffin [20]. 

We use a quantized-field representation with field op- 
erator '^(x), and canonical conjugate i^^(x): 



The Hamiltonian is 



(14) 



H= I dhW(xyh(x)1'(x) 
+ i I d'xd'yW(yyW(xyV(x,y)W(x)W(yX (15) 



where 



A(x) = -|^ V^+Vext(x)-)a. (16) 



where Ve^tix) is the external potential, and /jl is the 
chemical potential. We displace the field operator by its 
ground-state expectation: 

nx) = ilj(x)^<l>(xX (17) 

where ijj(x) is the displaced field operator, and 

<i>(x) = {nx)y (18) 

The function (t)(x) describes a nonuniform condensate, 
and ijj(x) annihilates particles not in the condensate. 

Variational trial states are taken to be the eigenstates 
of an effective Hamiltonian that is quadratic in the field 
operators, chosen to be a mean-field version of the ac- 
tual Hamiltonian. There are two variational functions 
p(xj) and A (x ,j ), which will turn out to be 
p(xj) = {ilj\x)ilj(y)), and zi(x j) = <^(x)^(y)). That 
is, p and A are, respectively, the direct and off-diagonal 
density correlation functions of the uncondensed parti- 
cles. 

The effective Hamiltonian is chosen to be 



H,,, = H'''^ + H''^ + H''^' + H'^\ 



(19) 



where //^^^ is independent of t/^, and H^^^ and //^^^ terms 
are respectively linear and quadratic in i^j: 

//^'^= \ d?x4>\x)h{x)(^{x) 
+ \^ d'xd'y4>'{y)4>\x)V{x^)<l>{x)4>iy\ (20) 



H^'^= I d^xip\x)h{x)4>{x) 

+ J d'xd'yV{x^)ilj\x) [c^(y)p(y,x) + 4>{x)p{y,y) 

+ c^*(y)D(x,j)], (21) 

H^^^= I d\ip\x)h{x)ip{x) 

+ J d^xd^jV(x J) [^^(y)^(x)/?(x J) 

+ lP\x)lP{x)R{y^)^ 

+ 1 Jd^xd^jV(xj) [^%)^^(y)D(xj) 

+ ^(x)^(y)D*(x,j)], (22) 

where 

R{x^) = p{x^) + 4>\x)4>{y\ 
D{x,y)^A{x,y) + 4>{x)<l>iy). (23) 

For arbitrary p and A , we can diagonalize the effec- 
tive Hamiltonian through a generalized Bogoliubov 
transformation: 



^{x) = ^\Un{x)7]n - VI{X)7]1^. 



(24) 



where j]n and j]l are annihilation and creation operators 
satisfying 



['J?«,'J?I] = 5m 



(25) 



To preserve the canonical commutators for ipix), we 
require 

E Vu,{x)u:{y) - v;(x)v„(y)] = 8\x - j). (26) 

n 

We now require //eff to have the form 



//eff = £'o + 2 ^nVlVn, 



(27) 



or, equivalently, 






(28) 
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These conditions determine the functions v^ , u^ and (p . 
The condition for cf) comes from the requirement that in 
//eff the terms hnear in ij/ix) vanish. The resuUing equa- 
tions are 

h{x)4){x) + J d?yV{x^) [(^(y)p(y^) + (^(x)p(y j) 

+ 4>'{y)A{x^) + \4>{yt4>{x)\ = 0, (29) 

€nUn(x) = h(x)ur^(x)-\- d^yV(x,y)[un(x)R(y,y) 

+ u,(y)R(y,x) - v„(j)D(x j)], (30) 

- €nVn(x) = h(x)vn(x) + d^yV(x ^) [Vn(x)R(y,y) 

+ Vn(y)R(y^) - u,(y)D\x^)l (31) 

To determine p and zl , let I ^) be the lowest eigenstate 
of //eff : 



H,f,\0) = Eo\0). 



(32) 



We require that the ground-state energy be at a mini- 
mum with respect to wriations in p and A : 



S{H) = 0, 



(33) 



where the brackets denote expectation value with re- 
spect to 1^). In calculating the expectation value above, 
we use Wick's theorem: 

{ilj\x)ilj'(y)ilj(x)ilj(y)) = {ilj\x)ilj(y)){ilj'(y)ilj(x)) 

+ mx)ilj(x)){ilj'(y)ilj(y)) + {ilj\x)ilj'(y)){ilj(x)ilj(y)). 

(34) 

The calculation is facilitated by noting that 5(//eff) = 0, 
since <P is an eigenstate. Thus we can use the equivalent 
condition 



8{H - //eff) = 0. 

The results are as follows: 



(35) 



p(x J) = {ilj\x)ilj(y)) = ^Vn(x)v:(yX 

n 

4(x J) = {il>{x)il>(y)) = - ^u„(x)v:(y). (36) 

n 

When these are substituted into Eqs. (29), (30), and 



(31), we have a system of coupled nonlinear integro-dif- 
ferential equations. 



5. The Uniform Case 

We put Vext = 0, assume V(x j ) = V(x — J ), and seek 
solutions with uniform density by putting 






(37) 



where Xk, y^ and (j)o can be taken to be real. The condi- 
tion [Eq. (26)] becomes xl — yl= \, which can be satis- 
fied by putting 



Xk = coshojt , 
yk=sinhak. 



(38) 



This parametrization simplifies the calculations. For ex- 
ample, the liquid-structure factor of the ground state can 
be represented in the form 



5jt = 1 + — (cosh 2ak — sinh 2aj^ — 1) 

n 



+ -T-jj 2 [sinh 2ap sinh 2a\k+p\ 

-\- sinh^ ap sinh^ oip+jti], (39) 

where Uq is the condensate density given through 

n = nQ-\- {2~^^ sinhVp . (40) 

p 

The Bogoliubov result is recovered by neglecting the/; 
sums, which is equivalent to putting Hq — n\ 



' V€, + 2nV,\ ' 



(41) 



where €k = ff'k^l2m, and V^ is the Fourier transform of 
the interaction potential. 

6. Reduction to 6 -Function Potential 

The equations in the last section are rather complex, 
and difficult to solve even numerically. To gain some 
insight, we reduce them to a simpler form by using a 
5-function potential. We shall show that this potential 
has no effect if self-consistency is enforced, at least in 
the case of a uniform system. We put 
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V(xj) = g8\x-yX (42) 

where 

g = 47T^^a/m. (43) 

The functions p and A now depend only on x : 

p(x) = {il/\x)ilj(x)) = ^Vn(x)vn(x) 

n 

A{x) = {il/(x)il/(x)) = - ^Un(x)v:(x). (44) 

n 

Note that p is the depletion of the condensate due to 
interactions. Equations (29), (30), and (31) reduce to 



[h + 2gp + g\<l>P]<l> + gA<l>^ = 0, 



(45) 



(h + 2gp + 2g(l)*4>)Un - g(A-\- (lp-)Vn = €Un, 

(h + 2gp+ 2g(^»v„ - g(^* + (A*')w« =- €Vn, (46) 

where we have suppressed the a: dependence of the func- 
tions. The dilute uniform hard-sphere gas is recovered 
by setting p = Zi = 0, for they give higher-order contri- 
butions. 

For p= A = 0, Eq. (45) is the "nonlinear Schrodinger 
equation" [21] [22], which has been widely used in 
numerical studies of the condensate and its excitations 
[23] . One should be aware, however, that p and A may be 
important, especially for the excitations. 

In the uniform case, with Vext = 0, we obtain, in the 
notation of the last section. 



tanh 2ak = 



g(n- p-\-A) 
€k-\-g(n- p- A)' 



(47) 



where €k = h^k^l2m, and the chemical potential }x has 
been eliminated in terms of the density n . The self-con- 
sistency conditions [Eq. (44)] then lead to 

1 r^ 

P = 4:^J ^kk^{fkV^k + g{n- p- A)\-\}, 
^="4^1 dkk%{n-p+A\ (48) 

where 

/, = [el + 2g€,{n -p-A)- Ag\n - p)AY"\ (49) 

The integrals are cut off at A, because otherwise A 
would be divergent. Solving for A from the second equa- 
tion, we obtain in the limit A ^^ 



The first equation then gives, in the limit A ^ oo, 

p = 0. (51) 

That is, there is no depletion of the condensate. It is then 
straightforward to show that the system is an ideal Bose 
gas. This result is reassuring; it shows that the method 
is able to verify that the 5 -function potential has no 
effect. 



7. Self-Consistent Field Method: Finite 
Temperatures 

To extend the variational approach to finite tempera- 
tures, we can use the eigenstates of ^eff to calculate the 
partition function, and then minimize the free energy 



F={H)- TS, 



(52) 



where the brackets now denote thermal average with 
respect to H^ff, and S is the entropy. De Gennes [19] 
argues that the quantity 



^eff - \^eff) ~ TS 



(53) 



is stationary, and implements the variational principle by 
minimizing F — Feff. But the assertion is correct only if 
S were defined with respect to the eigenvalue spectrum 
of ^eff- In fact, as we explain below, S is defined with 
respect to the energy levels (a\H\a), where la) are the 
eigenstates of ^eff- The finite- temperature results in Ref. 
[19] are therefore incorrect. 

A proper variational principle for the free energy [24] 
is based on the inequality 



2<ale-^''la)^Ee"^^""""^ 



(54) 



This shows that, to calculate the trial partition function, 
we must use the energy spectrum 



Wa = {a\H\a). 



(55) 



A - p-\-n = 0. 



(50) 



In contrast, the other scheme leads to a free Bose gas 
with single-particle energy 6„, the eigenvalues from 
Eqs. (30) and (31). Because of the simplicity, it is 
widely used in the literature; but it does not follow from 
a variational principle. 

As an example of the difference between the two 
schemes, the calculation of the partition of a hard- 
sphere Bose gas in Ref. [25] corresponds to using Eq. 
(55), and yields the virial coefficients in the gas phase. 
On the other hand, a calculation based on the other 
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scheme would give a free Bose gas above the transition 
temperature. 

8. Gaussian Variational Method 

We go to a representation in which the field operator 
^(x) and its conjugate have the forms 



nx) = 



1 

V2 



i^^^X): 



V2 



(p(x)-\- 
(p(x) - 



d(p{x) 

8 

8(p(x)_ 



(56) 



where (p(x) is a c -number function. The more familiar 
representation, in which ^ is diagonal, would be awk- 
ward to use in a non-relativistic case, in which ^ and 
i'^^ are canoncial conjugates. (A similar situation hap- 
pens in a relativistic fermion field obeying the Dirac 
equation.) 

The state of the system is described by a wave func- 
tional <P[(p]. We use a Gaussian form for a trial wave 
functional: 

^[cp] = C expj- i J d'xd'y[cp(x) 



(57) 



cpo(x)]G-\xj)[cp(y) - cpo(y)]l 



where C is a normalization constant and ^o and G(xj) 
are the variational parameters. Expectation values are 
given by functional integrals, for example 

{H) = j (Dcp)^\cp]H^[cp] (58) 



where ^ is normalized to unity. It is easy to show that 

(59) 



{nx)) = cpo(xx 

{<p(x)<p(y)) = G(xj>) -F (po(x)(po(y). 



The variational principle states that 






= 0, 



S<poix) 



= 0. 



(60) 



The equation obtained for G and <po will not be given 
directly. To make contact with the self-consistent field 
method, we quote the results 



p(jCJ')s<,/,(jc)V(y)) 



^G-\x^) + G(x^)-S\x-y) 



Aix^)^{iljix)ilj(y)) 
Ig-\x^)-G(xj) 



(61) 



For further analysis, it is convenient to introduce 
1 



X{XJ)-: 

Y(xj)^ 



"2V2 

1 

2V2 



G-^^^(xj) + 2G^^(xj) 



G-''\x^)-2G''^(xj) 



(62) 



it easy to show that p(x ^) = J d^zY(x ^)Y(zj>), and 
A(xj)=- J d^zX(x,z)Y(zj). Now expand Zand Fin 
terms of a basis [b^ix)}: 



X(xj) = J,XnK(x)b:(yX 

n 

Y{x^) = J,YJb„ix)b:{y), 

n 

where Z„ and F„ are real. We then find 

p{x^) = ^b^{x)b:{y)Yl 

n 

A{x,y)=-J,b,{x)b:{y)Xl 



(63) 



(64) 



Comparison with Eq. (36) leads to the identification 

(65) 



Un{x) = bn{x)Xn, 
Vn{x) = bn(x)Yn. 



This shows the equivalence between this method and the 
self-consistent field method. 

One of the advantages of the Gaussian method lies in 
the possibility of generalization to the "time- dependent 
Hartree-Fock" method [26,27], which has been used 
successfully in nuclear physics. One of us (P.T.) plans to 
discuss this in a separate publication. 
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